function PlotConcentrationDistribution_mut1_mut2(colonySet1, colonySet2, colonyBckgrnd)
%GETAUTOCORRELATION Summary of this function goes here
%   Detailed explanation goes here

%% Gather f and df data for suffix1
df1 = []; f1 = []; len1 = []; conc1 = []; totalCells1 = 0;
for j=1:length(colonySet1)
    cell = colonySet1{j}.cell;
    for i=1:length(cell)
        if strcmp(cell(i).fate, 'd')
            df1 = [df1 cell(i).flscnc.df];
            f1 = [f1 cell(i).flscnc.f];
            len1 = [len1 cell(i).size.l];
            conc1 = [conc1 cell(i).flscnc.fFilt];%./cell(i).size.lFilt];
            totalCells1 = 1 + totalCells1;
        end
    end
end

%% Gather f and df data for suffix2
df2 = []; f2 = []; len2 = []; conc2 = []; totalCells2 = 0;
for j=1:length(colonySet2)
    cell = colonySet2{j}.cell;
    f2col = [];
    for i=1:length(cell)
        if strcmp(cell(i).fate, 'd')
            df2 = [df2 cell(i).flscnc.df];
            f2 = [f2 cell(i).flscnc.f];
            len2 = [len2 cell(i).size.l];
            conc2 = [conc2 cell(i).flscnc.fFilt];%./cell(i).size.lFilt];
            totalCells2 = 1 + totalCells2;
        end
    end
end

%% Gather f and df data for suffixb (background)
dfb = []; fb = []; lenb = []; concb = []; totalCellsb = 0;
for j=1:length(colonyBckgrnd)    
    cell = colonyBckgrnd{j}.cell;
    for i=1:length(cell)
        if strcmp(cell(i).fate, 'd')
            dfb = [dfb cell(i).flscnc.df];
            fb = [fb cell(i).flscnc.f];
            lenb = [lenb cell(i).size.l];
            concb = [concb cell(i).flscnc.fFilt];%./cell(i).size.lFilt];            
            totalCellsb = 1 + totalCellsb;
        end
    end
end

%pdf
%.................................................
disp('Set1 statistics')
[f1 indx] = sort(f1);
len1 = len1(indx);
lb = f1(floor(0.17*length(f1)));
ub = f1(ceil(0.83*length(f1)));
averf1 = mean(f1)
sigma = std(f1) %(ub-lb)/2

noise = std(f1)/mean(f1)

disp('Set2 statistics')
[f2 indx] = sort(f2);
len2 = len2(indx);
lb = f2(floor(0.17*length(f2)));
ub = f2(ceil(0.83*length(f2)));
averf2 = mean(f2)
sigma = std(f2) %(ub-lb)/2

noise = std(f2)/mean(f2)

disp('SetBckgrnd statistics')
[fb indx] = sort(fb);
lenb = lenb(indx);
lb = fb(floor(0.17*length(fb)));
ub = fb(ceil(0.83*length(fb)));
aver = mean(fb)
sigma = std(fb) %(ub-lb)/2

disp('SetBckgrnd activity statistics')
dfb = sort(dfb);
lb = dfb(floor(0.17*length(dfb)));
ub = dfb(ceil(0.83*length(dfb)));
aver = mean(dfb)
sigma = std(dfb) %(ub-lb)/2

bckgrndLevel = aver + 3*sigma

disp(['Cells in set 1: ' num2str(totalCells1)]);
disp(['Cells in set 2: ' num2str(totalCells2)]);


%Obtain cdf for bursts
%................................................
range = -1000:100:7000;

df1 = df1(df1>=bckgrndLevel);
df2 = df2(df2>=bckgrndLevel);

cdf1 = zeros(1,length(range));
cdf2 = zeros(1,length(range));
norm1 = length(df1);
norm2 = length(df2);
for i=1:length(range), 
    cdf1(i) = sum(df1<=range(i))/norm1;
    cdf2(i) = sum(df2<=range(i))/norm2;
end


%% Concentration ratio
conc1 = sort(conc1);
conc2 = sort(conc2);
maximum = 20000;
res = 500;
ratio = zeros(1,maximum/res+1);
total1 = sum(conc1);
total2 = sum(conc2);
for i=maximum:-res:0
    thresh = i;
    indx1 = find(thresh <= conc1);
    indx2 = find(thresh <= conc2);
    if indx1
        indx = i/res + 1;
        ratio(indx) = (sum(conc2(indx2:end))/total2)/(sum(conc1(indx1:end))/total1);
    end
end

figure('Color', [1 1 1]);
plot(0:res:maximum, ratio, 'b-', 'Linewidth', 3)
line([-1 20000], [1 1], 'Color', 'k', 'LineStyle', ':');

m_concb = mean(concb);
line([1 1]*m_concb, [.01 1000], 'Color', 'k', 'LineStyle', '-.');
str1 = {'{\bf\langle}autofluorescence{\bf\rangle}'};
text(m_concb-500, 1.1, str1, 'FontSize', 13, 'Rotation', 90)

m_conc1 = mean(conc1);
line([1 1]*m_conc1, [.01 1000], 'Color', 'k', 'LineStyle', '-.');
str1 = {'{\bf\langle}mEGFP{\bf\rangle}'};
text(m_conc1-500, 1.55, str1, 'FontSize', 13, 'Rotation', 90)

m_conc2 = mean(conc2);
line([1 1]*m_conc2, [.01 1000], 'Color', 'k', 'LineStyle', '-.');
str1 = {'{\bf\langle}mGFPmut2{\bf\rangle}'};
text(m_conc2-500, 2.3, str1, 'FontSize', 13, 'Rotation', 90)

xlabel('\fontsize{18}Fluorescence Threshold, F_{0}')
ylabel('\fontsize{18}{\bfP}(F_{\fontsize{11}Mut2}>F_{0})/{\bfP}(F_{\fontsize{11}EGFP}>F_{0})')

xlim([0 15000])
ylim([0.85 10])

set(gca, 'YScale', 'log');
set(gca, 'FontSize', 18);


%% .........................................................
df1 = df1(df1>-1000);
df2 = df2(df2>-1000);
dfb = dfb(dfb>-1000);
median(df1)
median(df2)
median(dfb)

[f1,x_f1] = ksdensity(f1, 'width', 200, 'npoints', 400);
[f2,x_f2] = ksdensity(f2, 'width', 200, 'npoints', 400);
[fb,x_fb] = ksdensity(fb, 'width', 200, 'npoints', 400);

figure('Color', [1 1 1]);
hold on;
plot(x_f1*.0001, f1, '-', 'LineWidth', 6, 'Color', [16 147 11]/255)
plot(x_f2*.0001, f2, '-', 'LineWidth', 6, 'Color', [72 244 66]/255)

xlim([0.1 2]);
ylim([0 0.00025]);

box on;

set(gca,'box','off')
isholdonque = ishold; 
hold on 
ax = axis; 
plot(ax(2)*[1,1],ax(3:4),'k','linewidth',0.5) 
plot(ax(1:2),ax(4)*[1,1],'k','linewidth',0.5) 
if isholdonque == 0 
hold off 
end

xlabel('\fontsize{19}Fluorescence Signal per cell (au)')
ylabel('\fontsize{20}Frequency')

averf1 = averf1*.0001;
averf2 = averf2*.0001;
set(gca, 'FontSize', 15);

hleg1 = legend('Slow FP','Fast FP');
set(hleg1,'FontSize',20);
legend(hleg1, 'boxoff')


end
